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ABSTRACT 


The  mechanical-graphical  “peel-off*  method  and  Marquardt’s  composite  Gauss- 
Newton  and  gradient  iterative  method  were  programmed  for  the  Philco  2000,  a  16K 
asynchronous  digital  computer.  Both  programs  were  coded  in  the  Philco  Algebraic 
Programming  Language  (ALTAC)  using  single-precision  floating-point  arithmetic. 

Background  material,  flow  charts,  flow  chart  descriptions,  subprogram  usage, 
computer  memory  requirements,  and  illustrative  numeric  examples  of  the  analyses 
of  both  simulated  and  empirical  data  are  given.  Each  sample  of  simulated  data 
possessed  an  error  component;  the  effects  of  an  asymptote,  in  several  instances,  were 
included  during  the  generation  of  the  data.  Dog  lung  nitrogen  washout  activity 
experiments  were  the  source  of  the  empirical  data. 


FITTING  MULTI  COMPONENT  EXPONENTIAL  DECAY  CURVES  BY  DI6ITAL  COMPUTER 


I.  INTRODUCTION 

The  principal  purpose  of  this  paper  is  to  discuss  the  automation  of  two  non¬ 
linear  parameter  estimation  procedures — the  classical  “peel-off”  and  Marquardt’s 
algorithm  (2).  As  a  mechanical-graphical  method,  the  “peel-off”  procedure 
yields  good  results  for  parameter  estimation  in  a  hypothesized  mathematical 
model  of  a  linear  combination  of  exponential  functions.  But  the  length  of  time 
required  in  the  application  of  this  method  is,  in  general,  too  long.  Automation 
of  the  “peel-off’  method  rectified  the  time  problem,  but  the  computer-produced 
parameter  estimates  turn  out  to  be  somewhat  inferior  in  accuracy.  A  parameter- 
estimate-refining  program  was  written  to  improve  the  “peel-off”  estimates.  The 
refining  process  was  accomplished  through  the  adaptation  of  an  algorithm  that 
was  described  by  Marquardt  (2)  for  obtaining  least-squares  estimates  of  non¬ 
linear  parameters.  Two  iterative  methods,  classical  Gauss-Newton  and  gradient, 
were  combined.  This  combination  yielded  an  iterative  method  with  strong  con¬ 
vergence  properties  and  a  compromise  between  two  levels  of  rapidity  of  con¬ 
vergence.  Both  curve-fitting  programs  were  coded  in  the  Philco  Algebraic 
Programming  Language  (ALTAC)  (3),  using  single-precision  floating-point 
arithmetic.  The  computational  mode  has  a  range  from  slightly  more  than  iO"00 
to  slightly  less  than  IO-*1"  and  an  accuracy  of  ten  significant  digits.  The  pro¬ 
grams  were  tested  on  simulated  as  well  as  on  dog  lung  nitrogen  washout  data. 

Four  sections  follow  this  introductory  section.  The  mathematical  model  used 
in  fitting  a  linear  combination  of  exponentials  is  briefly  discussed  in  section  II. 
The  reader  is  referred  to  Danford  (1)  for  a  genera)  discussion  of  exponential 
model  equations.  Background  material,  flow  chart,  flow  chart  description,  sub¬ 
program  usage,  computer  memory  requirements,  and  numerical  examples  of  the 
analyses  of  simulated  data  (generated  with  and  without  the  effects  of  an  asymp¬ 
tote  and  with  error)  and  dog  lung  nitrogen  washout  data  pertinent  to  the  “peel- 
off  method  are  covered  in  section  III.  Material  similar  to  that  of  section  III,  on 
the  composite  Gauss-Newton  and  gradient  method,  is  to  be  found  in  section  IV. 
Section  V  is  concerned  with  comments  on  such  selected  material  as:  adapting 
the  computer  programs  to  meet  the  user's  requirements;  possible  places  for  im¬ 
provement  in  the  “peel-off'  method  program;  the  consequence  of  failing  to  use 
the  best  possible  estimate  of  the  asymptote;  the  use  of  smoothed  data  to  improve 
parameter  estimation;  the  importance  of  refining  the  preliminary  parameter 
‘  estimates  produced  by  the  “peel-off'  method  program ;  and  the  area  in  which  the 
use  of  double-precision  floating-point  arithmetic  may  become  necessary  in  the 
Gauss-Newton  and  gradient  iterative  process. 


II.  MATHEMATICAL  MODEL 


Our  main  results  concern  the  estimation  of  the  parameters  in  the  model 

Y(x)  =  y(x)  +  «(x)  (1) 

=  «0  +  I«„  exp(-PDx)  +  t(x) 

m— 1 

N 

=  »o  +  2  +  •(*)  » 

m-«I 

using  the  data  points  (Xi ,  Y,),  i  =  1,2 . L,  where 

y(x)  :  True  value  of  Y  at  x, 

»0  :  Constant  term  or  asymptote, 

•m,  pm  :  Model  parameters  >  0,  m  =  1,2 , . . . ,  N, 

«(x)  :  Error  term:  NID  (0,v)  for  each  x  and 

•  •  t  y(x), 

p  :  Positive  number;  100p  can  be  viewed  as  a  percent  error, 

L  :  Total  number  of  data  points, 

N  :  Number  of  exponential  components. 


III.  “PEEL-OFF”  METHOD 


Meet  anical-graphical  version 

The  “peel-off  method  for  parameter  estimation  has  been  used  for  a  con¬ 
siderable  period  of  time  without  undergoing  any  major  modification.  Its  users 
during  this  time  have  probably  numbered  in  the  hundreds  and  this  extensive 
usage  has  no  doubt  been  due  to  the  absence  of  a  superior  technic.  Application 
of  the  method  is  quite  easy.  Simple  tools  such  as  a  pencil,  semilogarithmically 
scaled  paper,  and  five-place  common  logarithm  tables  are  needed  plus  “good” 
judgment  on  the  part  of  the  analyst.  The  method  can  become  very  tedious 
when  many  parameter  estimates  per  experiment  are  required  or  there  are  many 
experiments  to  be  analysed.  One  must  not  fail  to  mention  disadvantages  of  a 
graver  nature:  personal  bias  may  be  present  that  will  affect  the  fit;  the  use 
of  imperfect  semilogarithmic  paper  may  contribute  to  the  bias  in  the  estimates; 
the  errors  of  estimation  may  be  cumulative  in  nature;  and  finally  there  is  the 
problem  of  securing  an  estimate  of  error. 

A  simplified  summary  of  the  sequence  of  steps  taken  by  the  data  analyst  in 
applying  the  “peel-off  method  is  ss  follows: 

Step  1.  Obtain  an  estimate  of  the  asymptote  a«  by  any  productive  means 
available. 


Step  2.  Remove  the  effects  of  the  asymptote  from  the  data  by  decreasing 
each  datum  by  the  magnitude  of  the  asymptote. 
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Step  3.  Plot  the  residuals,  obtained  in  step  2,  versus  x,  on  h-cycle  base  10 
semilogarithmic  paper. 

Step  4.  Fit  a  straight  line  by  “eye”  to  as  many  points  as  are  judged  a 
“good”  fit.  Start  fitting  the  points  in  the  right-hand  end  of  the  plotted  curve. 
Furthermore,  after  the  straight  line  is  fitted,  extrapolate  back  to  the  semi- 
logarithmic  axis  to  produce  an  estimate  of  parameter  a»,  ;  the  estimate  of  the 
associated  pm  is  obtained  through  the  use  of  a  simple  modified  form  of  the 
analytic  expression  for  the  fitted  line. 

Step  5.  Terminate  the  "peeling”  process  if  there  are  fewer  than  four  points 
to  fit ;  otherwise,  proceed  to  step  6. 

Step  6.  Obtain  a  new  set  of  residuals  by  subtracting  the  effects  of  the  com¬ 
ponent  fitted  in  step  4,  using  original  units,  from  the  points  that  were  not  in¬ 
cluded  in  the  fit. 

Step  7.  Plot  the  residuals,  obtained  in  step  6,  versus  x,  on  semilogarithmic 
paper.  Return  to  step  4. 

Figure  1  depicts  an  example  of  component  “peeling”  by  means  of  the  seven 
aforementioned  steps.  In  this  particular  instance,  four  components  were  ex¬ 
tracted  from  the  data  after  the  effects  of  a  constant  (asymptote,  estimated  for 
this  case  as  two-thirds  of  the  last  recorded  Y  value)  had  been  removed.  After 
the  fit  of  each  respective  straight  line,  the  parameter  estimate  of  a  is  obtained 
by  taking  */«„.,  (since  the  original  Y  values  were  multiplied  by  100  before  plot¬ 
ting)  the  value  read  off  the  semilogarithmic  axis.  The  simple  formula 

ft  —  —  (log?  —  loga)lnl0/x  (X) 

was  then  used  to  compute  the  parameter  estimate  of  p  when  the  associated 
estimate  of  a  became  available.  The  symbols  log  and  In  stand  for  common  and 
natural  logarithms  respectively ;  the  single  point  (J,1?)  used  is  any  point  on  the 
fitted  straight  line.  The  estimate  of  the  parameter  »  is  exp(-jj).  Observe  the 
estimates,  by  component,  in  the  legend  of  the  graph. 

Digital  computer  version 

We  now  proceed  to  describe  the  programmed  version  of  the  “peel-off*  meth¬ 
od.  Our  description  will  cover  input,  general  computational  steps,  fitted  straight- 
line  acceptance  procedure,  additional  computations,  and  output.  Following  this 
there  will  appear  a  flow  chart  of  the  program  along  with  its  description,  sub¬ 
program  usage,  computer  memory  requirements,  and  illustrative  numeric  ex¬ 
amples  of  program  yield  in  the  analysis  of  simulated  and  empirical  data. 

Input: 


1.  Program  parameters  (from  a  punched  card). 

a.  EXP:  Experiment  number. 

b.  SAMP:  Sample  number. 
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c.  L:  Sample  size. 

d.  N :  Number  of  components  expected  to  be  fitted. 

e.  a0:  Estimate  of  the  asymptote. 

2.  Data  (from  punched  cards) . 
x,Y  values. 

Remarks :  When  the  asymptote  is  absent  from  the  mathematical  model,  its 
estimate  is  zero ;  otherwise,  obtain  a  non-zero  estimate  that  is  based  on  the  ex¬ 
perimenter’s  personal  experience,  a  graphical  approach,  or  any  other  means  that 
will  provide  a  satisfactory  preliminary  estimate  of ‘this  parameter. 

Genera]  computational  steps: 

Step  1.  Obtain  an  estimate  of  the  asymptote. 


Step  2.  Remove  the  effects  of  the  asymptote  from  the  data. 

RY,  =  Y,  —  ,  1  =  1,2 . L.  (8) 

Step  3.  Rearrange  the  residuals  obtained  in  step  2, 

(*L+1-,.RYl  +  1_l)-*(Z,Jt1)  (4) 

and  transform  the  residuals  R,  logarithmically, 

G,  =  Ln(R,),  i  =  12 . L.  (5) 

Step  4.  Fit  a  straight  line,  by  principle  of  least  squares,  to  K  points  (Z( ,  G,). 
Start  fitting  with  i  =  1.  (See  line-fitting  acceptance  procedure  following  step  7 
below.) 

Use 

G  —  Ln(ft)  —  yZ  ,  Z|  —  Z  —  Zg  (6) 

to  obtain  estimates  of  a,  fi  for  the  mth  component. 

Step  5.  Terminate  the  “peeling”  process  if  there  are  fewer  than  four  points 
remaining  after  deleting  the  K  points  used  in  step  4 ;  i.e.,  we  must  have 

NR  =  (NPR  —  K)  <  4  U> 


where 

NPR:  Total  cumber  of  poinU  available  for  fitting  proceee  in  step  4  (initially  NPR  =  L). 

Otherwise  proceed  to  step  6. 

Step  6.  Obtain  a  new  set  of  residuals  by  removing  the  effects  of  the  mth 
component  “peeled”  in  step  4. 

RPj  =  Rj  -  s,  «P( -  ?«Z,>  .  j  =  R  +  ijt +  NPR  .  (t) 
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Step  7  Delete  the  K  elements  used  in  step  4  by  repositioning  the  G,R,Z- 
arrays 


ZJ  -  ZK+J 

(9) 

Rj  =  RPr+j 

(10) 

Gj  =  Ln(Rj)  , 

j  =  1,2 . NR  . 

(ID 

Set  NPR  =  NR  and  return  to  step  4. 

In  the  acceptance  of  the  existence  of  a  “real”  component  which  will  be  dis¬ 
cussed  next,  the  motivation  for  all  the  steps  will  not  be  given,  for  the  whole 
procedure  is  subjective.  The  automation  of  this  scheme  is  based,  primarily,  on 
our  experience  with  one  kind  of  data — nitrogen  washout  of  the  lung.  It  is 
hoped,  however,  tl.at  the  procedure  will  be  fairly  general  and  that  the  acceptance 
tests  to  be  described  below  will  have  some  face  validity  to  the  reader.  A  slight 
note  is  occasionally  made  to  motivate  a  test,  but  there  is  no  claim  for  this  being 
the  best  procedure  that  can  be  developed  for  Ihe  “peel-off”  method.  It  has 
undergone  several  revisions  in  our  laboratory  and  must  be  considered  as  a  best 
effort  at  this  point  in  time. 


Acceptance  procedure: 

A  fitted  least-squares  line,  which  implies  the  existence  of  a  real  component, 
except  for  the  “last  component,”  is  accepted  if  the  following  tests  are  satisfied 
in  the  stated  order: 

1.  Runs  test. 

The  differences  for  the  next  four  successive  points  are  all  positive. 

DlFFj  =  Rj  -«m  exp(-^mZJ) 

=  Rj  -  YEST,  >  0,  j  =  K  +  1,K  +  2 . K  +  4  .  (12) 


2.  Beta  test 

The  slope  of  the  least  squares  line  fitted  to  the  (m-l)st  interval  is  less 
than  the  slope  of  the  line  fitted  to  the  mth  interval. 

Pm  >  Pm-l  .  (13) 


where 

K  -  0.0 


3.  F  test. 

The  ratio  is 

DEVINC/4 

F  - - >  TOL  , 

DEVFIT/(K  —  2) 


(14) 
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where 

K  +  4 

DEVINC  =  2  (R.  -  YEST,)2 

j  =  K  +  1 

K 

DEVPIT  =  2  (Rj  -  YEST,)2 
j  =  1 

and 


K:  Number  of  points  in  the  interval  of  fit, 

TOL:  Upper  1%  point  of  the  Snedecor-Fisher  (F)  distribution  with 
4  and  K-2  degrees  of  freedom,  for  K  ^  32, 

:  4,  for  K  >  32  . 


4.  Remaining  res.duals  test. 

The  differences  for  all  successive  residuals  not  included  in  the  fit  are  posi¬ 
tive. 

RES,  =  R,  -  YEST,  >0,  j  =  K  +  1,K  +  2 . NPR  ,  (16) 

where 


NPR:  Number  of  residuals  yet  to  be  fitted. 


Note  that  test  1  is,  in  general,  a  subtest  of  this  test.  It  is  used  as  a 
quick  or  preliminary  test.  This  more  extensive  test  4  is  needed  to  overcome  cer¬ 
tain  abnormal  deviations  which  test  1  will  not  detect. 


6.  Alpha  test. 

The  estimate  of  am  +  ,,  computed  in  the  fitting  of  the  next  four  successive 
noints,  must  be  less  than  or  equal  to  2Y,,  i.e., 

Vi~2Yi  •  (1«) 

Here  Y,  is  the  first  observation  of  the  original  data. 


This  test  is  not  as  obvious  as  the  others.  It  has  been  empirically  deter¬ 
mined  to  guard  against  including  too  many  points  in  the  interval  of  fit  of  the 
component  that  is  being  currently  fitted. 


The  logic  of  the  fitting  process  is  such  that  from  one  up  to  and  including 
a  predetermined  number  of  components  is  fitted.  Furthermore,  four  points  is 
the  initial  number  as  well  as  the  least  number  of  points  which  is  included  in  a 
fitted  component.  Also  for  greater  clarity,  two  distinct  categories  of  fittings 
need  to  be  considered :  fitting  the  first  through  the  next  to  the  last  component 
and  fitting  the  last  component.  Here  "last  component”  is  defined  as  a  state 
that  exists  when  either  one  of  the  following  two  conditions  holds; 
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1.  The  predetermined  number  N  of  components  expected  to  be  fitted  has 
been  reached. 

2.  NPR  -  K  ^  4. 


For  either  condition  1  or  2  holding,  all  acceptance  procedure  tests  are  ignored 
when  we  fit  the  last  component. 


1.  Fitting  the  first  through  the  next  to  the  last  component: 


If  test  1  is  not  satisfied:  Then  either  there  was  or  was  not  a  run  of  four 
minus  signs.  In  the  case  of  four  minus  signs,  shift  the  G,R,  and  Z-arrays  such 

that  for  each  array  the  (i  -f-  l)st  element  replaces  the  ith  element:  refit  a  '■ 

least-squares  line  to  interval  of  K  points  of  repositioned  G,  Z-arrays.  Repeat  this  >. 

point  deletion  action  (referred  to  as  “creeping”)  and  refitting  process  until  | 

test  1  is  satisfied.  In  the  other  case,  an  additional  point  is  included  in  the  .in-  i 

terval  and  the  least-squares  line  is  refitted.  The  process  is  continued  until  test  1  I1 

is  satisfied.  This  test  is  admittedly  arbitrary  and  can  be  altered  so  that  one  is 
more  or  less  certain  of  detecting  the  beginning  of  another  component. 

If  test  2  is  not  satisfied:  Then  “creep”  in  the  manner  described  above  » 

for  test  1.  The  “creeping”  and  line  refitting  process  continues  until  test  2  is  | 

satisfied.  j 

If  test  3  is  not  satisfied :  Proceed  in  the  same  manner  as  for  test  1  for  S 

the  case  in  which  a  run  of  four  minus  signs  did  not  occur.  j 


If  test  4  is  not  satisfied :  Repeat  action  similar  to  that  taken  for  test  8.  i 

i 

\ 

If  test  5  is  not  satisfied:  This  indicates  that  too  many  points  were  t 

used  in  the  fitting  of  the  least-squares  straight  line.  Proceed  to  reduce  the  ^ 

number  of  points  in  the  interval  of  fit  one  point  at  a  time  and  use  the  following  » 1 

criteria  for  terminating  the  process:  , 

I 


Test  1,  2,  or  4  fails — restore  one  point  to  the  reduced  interval;  refit  the 
new  interval;  accept  the  fitted  line, 


Test  1,  2,  4  and  5  hold — accept  the  line  fitted  to  the  reduced  interval, 


Number  of  points  in  the  interval  has  been  reduced  to  only  four— accept 
the  line  fitted  to  these  four  points. 


< 


; 

I 
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2.  Fitting  the  last  component : 


The  program  will  fit  the  last  four  points  in  the  sequence  of  residuals 
that  remains.  Again,  fitting  exactly  four  points  is  arbitrary,  but  it  seems  to 
give  one  essentially  unbiased  estimates. 


Additional  computations : 

1.  Compute  a  mean  square  ratio  for  each  component 


SSRn,  = 


Km 

2  [<Rj  -  YESTj)/YESTj]2 


K„,  -  2 


»  m  —  1,2,  * . . ,  M  , 


where 


K,„:  Number  of  points  used  in  fitting  the  mth  component, 
M  :  Total  number  of  components  fitted. 


(17) 


2.  Compute  omegas 

=  exp(  -£„,),  m  =  1,2, . . . ,  M  . 

3.  Compute  estimated  Y  values 

M 

YE,  =  £0  +  2  amexp(  -  ?B1x,)  ,  i  =  1,2 . L  . 


4.  Compute  ratios 

RA,  =  (Y,  -  YE,)/YE,  ,  i  =  1,2 . L  . 

5.  Compute  cumulative  sums  of  squared  ratios 

CRATIO,  =  2  RA*.  ,  i  =  u . L. 

J-i 


6.  Compute  overall  “unrefined”  mean  square  ratio 


PM  SR  = 


CRATIO,. 
L  -  NP" 


(18) 


(19) 


(20) 


(21) 


(22) 
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where 


NP  = 


2M  for  asymptote  not  estimated, 
2M  +1  otherwise. 


Output: 

1.  Print — 

a.  Experiment  number. 

b.  Sample  number. 

c.  For  each  data  point: 

(1)  x,  Y  value. 

(2)  Estimated  Y  value. 

(3)  Ratio. 

(4)  Cumulative  sum  of  squared  ratios. 

d.  Model:  Model  1  for  no  asymptote. 

Model  2  otherwise. 

e.  Sample  size. 

f.  Number  of  components  expected  to  be  fitted. 

g.  Number  of  components  actually  fitted. 

h.  Overall  “unrefined"  mean  square  ratio. 

i.  Estimate  of  constant. 

j.  For  each  component  fitted : 

(1)  Range  for  x  values  included  in  the  fit. 

(2)  Estimate  of  alpha,  beta,  and  omega. 

(3)  Mean  square  ratio. 


2.  Write  on  output  tape— 


a.  Experiment  number. 

b.  Sample  number. 

c.  Model. 

d.  Sample  size. 

e.  Number  of  parameters. 

f.  Overall  “unrefined"  mean  square  ratio. 

g.  The  constant  and  an  alpha,  beta,  and  omega  for  each  component  fitted. 

h.  Data  points. 


FLOW  CHART  I 

FLOW  CHART  OF  THE  "PEEL-OFF1'  METHOD 
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IS 


Description  of  flow  ehsrt  I 


Box  1*4:  Self-explanatory. 

5:  The  following  arrays  are  cleared: 
a:  Parameter  array. 

0:  Parameter  array. 

«:  Parameter  array. 

SSR:  Mean  squared  ratios  for  the  components. 

FL:  Left-end  points  of  the  intervals  of  component  fit. 

FR:  Right-end  points  of  the  intervals  of  component  fit. 

X:  Abscissae  of  the  data  points. 

Y :  Ordinates  of  the  data  points. 

G:  Natural  logarithmic  transformation  of  R-array. 

R:  Residuals  stored  in  reverse  order. 

Z:  X-array  stored  in  reverse  order. 

<:  Self-explanatory. 

7:  Residuals  obtained  by  removing  the  effects  of  an  asymptote,  making  natural 
logarithmic  transformation  of  residuals,  and  developing  G,  R,  and  Z-arrays. 

8:  Initialize  component  number  counter  m  and  number  of  points  remaining  to  be 
fitted — gage  word  NPR. 

9:  Increment  the  component  number  counter  in,  initialise  the  number  of  points  to 
be  fitted  (K  =  4),  and  set  overfit  switch  (SWO)  off. 

19:  If  the  component  is  the  last  component  to  be  fitted,  go  to  LAST  ONE — box  12. 

11 :  If  a  sufficient  number  of  points  remain  for  a  normal  fit,  go  to  FIT— box  14. 

12:  Set  up  G,  R,  and  Z-arrays  so  that  the  last  four  points  can  be  fitted;  set  K  =  4  and 
NPR  =  4. 

19:  Fit  last  component;  go  to  ACCEPT — box  88. 

14:  A  least-squares  straight  line  is  fitted  to  K  points. 

IS:  Compute  differences  for  points  K  -f  1  to  X  +  4. 

II:  If  all  differences  are  negative,  go  to  CREEP— box  19. 

17:  If  all  differences  are  not  positive,  go  to  NO  RUN— box  89. 

18:  Beta  test  If  test  holds,  go  to  box  81. 

19:  Decrement  the  number  of  points  remaining,  NPR. 

91:  Shift  the  G,  R,  Z-arrays  such  that  for  each  array  the  (i  +  l)st  element  replacee 
the  ttfc  element;  go  to  PTS  LEFT— box  11. 

21:  If  refittuig,  go  to  CK  RES — box  28. 
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22-27:  F-test.  If  test  fails,  go  to  NO  RUN — box  29. 

28:  If  all  points  not  included  in  the  fit  lie  above  the  fitted  line,  go  to  CK  ALPHA — 
box  30. 

29:  Increment  counter  k;  go  to  PTS  LEFT — box  11. 

30-32:  Alpha  test.  If  test  holds,  go  to  ACCEPT — box  M. 

33:  If  the  overfit  switch  (SWO)  is  on,  go  to  box  86. 

34:  Overfit  switch  (SWO)  was  off — set  it  on;  go  to  TEST  K — box  30. 

35:  If  the  reduced  interval  of  fit  is  rejected,  accept  currently  fitted  component; 
go  to  ACCEPT — box  38. 

30:  If  it  is  impossible  to  reduce  the  current  interval  of  fit  by  one  point,  go  to 
ACCEPT— box  38. 

37:  Current  interval  of  fit  can  be  reduced — reduce  counter  k  by  one  and  set  KP 
equal  to  new  value  of  counter  k;  go  to  FIT — box  14. 

38:  Compute  sum  of  squared  ratios  for  the  K  points  fitted.  Save  various  quantities 
associated  with  the  mth  fitted  component  and  compute  the  number  of  points 
remaining. 

39:  If  no  points  remain,  go  to  OUTPUT — box  41. 

40:  Compute  residuals,  store  in  R-array;  compute  natural  logarithms  of  residuals, 
store  in  G -array;  reposition  Z-array;  go  to  NEW  COMP — box  9. 

41-43:  Compute  the  number  of  parameters  fitted  and  set  np  model  number. 

44-46:  Print  output: 

1.  EXP.  SAMP. 

2.  X,  Y,  YE,  RA,  CRATIO,  for  each  data  point. 

3.  EXP,  SAMP.  MODEL,  L,  N,  M,  PMSR. 

4.  «0;  FL,  FR,  SSR, «,  (1, 6  for  each  component  fitted. 

40:  Write  output  tape: 

1.  EXP,  SAMP,  MODEL,  L,  M,  NP,  PMSR. 

2.  ft,  and  estimates  t,  0,  S  for  eact  omponent  fitted. 

3.  X,  Y  values. 

Go  to  READ — box  3. 

Subprogram  usage 

Function  subprogram*  (the  first  two  are  standard  library  functions)  listed 
below  proved  helpful: 

1.  EXPF.  Argument:  A  (mention  of  expression  A).  Function:  computes 
the  value  exp  (A). 
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2.  LOGF.  Argument :  A  (location  of  expression  A).  Function:  computes 
the  value  In  (A). 

3.  YEST.  Argument:  I  (location  of  subscript  I).  Function:  computes  the 
ith  estimated  Y  value  using  the  estimates  fim  and 

A  subroutine  subprogram  needs  to  be  mentioned : 

FITLINE.  Argument:  A  (location  for  storing  a  ),  B  (location  for  stor¬ 
ing  ti  ),  C  (location  of  first  element  in  Z-array),  D  (location  of  first  element  in 
G-array),  and  J  (location  of  the  number  of  points  in  the  fit).  Function:  com¬ 
putes  estimates  of  the  a  and  /9  by  fitting  a  least-squares  line  to  the  natural 
logarithms  of  J  points. 


Memory  requirements 
Program 

A 

do 

A 

% 

A 

te>l 

X-array 
Y-array 
G-array 
R-array 
Z-array 
Range  limits 
Component  MSR’s 
Other 


About  2,680  words 
1  word 
M  words 
M  words 
M  words 
L  words 
L  words 
L  words 
L  words 
L  words 
2M  words 
M  words 
150  words 


Total :  2,831  -f  5L  -f  6M  words  approximately 


Examples 

Simulated  and  empirical  data  served  as  input  to  the  “peel-off”  method  com¬ 
puter  program.  The  generation  of  the  simulated  data  was  accomplished  through 
the  use  of  the  second  equation  of  1;  the  error  term  «(; >.  for  each  value  of  x, 
was  produced  by  means  of  a  subroutine  that  generated  pseudo-random  normal 
deviates  with  mean  —  0,  standard  deviation  -  #>y(x) ;  the  generator  of  the 
pseudo-random  numbers  used  in  the  normal  deviate  generation  was  of  the 
multiplicative  congruential  type.  The  empirical  data  was  restricted  to  dog  lung 
nitrogen  washout  data.  Program  yield  for  the  cases  considered  is  presented  in 
tables  I  to  ill  for  simulated  data  analysis  and  tables  IV  to  VI  for  empirical  data 
aiialysis. 
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TABLE  I 


Preliminary  parameter  estimates 
(Simulated  data) 


Sample  size:  SO 
Value  of  p :  0.0005 

Square  root  of  overall  MS  ratio:  0.0072 

True  value  of  o„ :  0.005  Estimate  of  a(l:  0.006* 


Component 

No. 

MS  ratio 

Range  of  x  values 
in  interval  fitted 

a 

P 

O) 

1 

12920  x  10  -10 

30.0  -  38.0 

0.020 

0.036 

0.965 

True 

Value 

0.020 

0.030 

0.970 

2 

296  X  10-i« 

1.0  -  4.0 

0.698 

0.367 

0.700 

True 

Value 

0.700 

0.367 

0.700 

•Estimate  of  a„  arbitrarily  chosen  in  the  neighborhood  of  the  true  value. 


TABLE  II 

Preliminary  parameter  estimates 
( Simulated  data ) 


Sample  size:  76 
Value  of  p:  0.0010 

Square  root  of  overall  MS  ratio:  0.0433 

True  value  of  None  Estimate  of  a„:  None 


Component 

No. 

MS  ratio 

Range  of  x  values 
in  interval  fitted 

a 

P 

(It 

1 

856  x  10-i» 

71.0  -  76.0 

0.065 

0.036 

0.965 

True 

Value 

0.040 

0.030 

0.970 

2 

31526  x  10- 111 

13.0  -  21.0 

0.398 

0.103 

0.902 

True 

Value 

0,400 

0.094 

0.910 

3 

109720  x  10-"' 

1.0  -  4.0 

0.180 

0.481 

0.618 

True 

Value 

0.200 

0.431 

0.650 

TABLE  III 


t  Preliminary  parameter  estimates 
(Simulated  data) 

Sample  size:  42 
Value  of  p:  0.0030 

Square  root  of  overall  MS  ratio:  0.0125 

True  value  of  oa:  0.050  Estimate  of  a0:  0.060* 


Component 

No. 

MS  ratio 

Range  of  x  values 
in  interval  fitted 

a 

0 

U) 

1 

1804  x  10-B 

32.0  -  40.0 

0.024 

0.084 

0.919 

True 

Value 

0.050 

0.105 

0.900 

2 

1968  X  lO-8 

18.0  -  28.0 

0.042 

0.162 

0.850 

True 

Value 

0.100 

0.357 

0.700 

3 

39235  X  10-s 

13.0  -  16.0 

0.006 

0.205 

0.814 

True 

Value 

0.200 

0.693 

0.600 

4 

430  x  10-« 

1.0  -  4.0 

0.631 

0.720 

0.487 

True 

Value 

0.400 

1.204 

0.300 

•Estimate  of  a0  simply  chosen  equal  to  the  true  value. 


TABLE  IV 

Preliminary  parameter  estimates 
(Empirical  data) 

Sample  size:  29 

Square  root  of  overall  MS  ratio:  0.0876 
Estimate  of  a0:  0.005* 


Component  ! 
No. 

MS  ratio 

Range  of  x  values 
in  interval  fitted 

a 

0 

« 

1 

2018  x  10-» 

26.0  -  29.0 

0.012 

0.080 

0.923 

2 

25218  x  10-» 

1.0  -  4.0 

0.890 

0.358 

0.699 

•Estimate  of  alt  based  on  a  visual  inspection  of  plotted  curve. 


TABLE  V 

Preliminary  parameter  estimates 
(Empirical  data) 


Sample  size:  49 

Square  root  of  overall  MS  ratio:  0.0665 
Estimate  of  a0:  None 


Component 

No. 

MS  ratio 

Range  of  x  values 
in  interval  fitted 

« 

0 

W 

1 

186  x  10- 

46.0  -  49.0 

0.016 

0.014 

0.986 

8 

4586  X  10-« 

87.0  -  43.0 

0.172 

0.168 

0.860 

8 

1678  x  10"1 

i 

1.0  -  4.0 

0.617 

0.362 

0.696 

TABLE  VI 

Preliminary  parameter  estimates 
(Empirical  data) 


Shmple  size:  140 

Square  root  of  overall  MS  ratio:  0.5504 
Estimate  of  a0:  0.010* 


Component 

No. 

MS  ratio 

Ranges  of  x  values 
in  interval  fitted 

a 

P 

Cl) 

1 

268  X  10-1° 

105.0  -  117.0 

0.006 

0.005 

0.995 

2 

800  X  10-6 

77.0  -  88.0 

0.013 

0.028 

0.972 

3 

2705  X  10-6 

1.0  -  4.0 

0.621 

0.260 

0.771 

•Estimate  of  a„  based  on  a  visual  inspection  of  plotted  curve. 


IV.  COMPOSITE  GAUSS-NEWTON  AND  GRADIENT  METHOD 
Marquardt’s  algorithm 

A  composite  Gauss-Newton  and  gradient  method  has  been  programmed  to  re¬ 
fine  the  preliminary  estimates  of  the  parameters  in  our  mathematical  model  that 
the  “peel-off”  program  provides.  It  is  an  itt  rative  method  that  was  designed  to 
eliminate  the  inadequacies  of  straight  Taylor-series  methods  and  gradient  meth¬ 
ods  but  at  the  same  time  retain  the  “good”  properties  of  both  approaches.  Ex¬ 
amples  of  good  properties  are:  convergence  greatly  accelerated  when  close 
proximity  to  converged  values  is  attained  and  region  of  convergence  broader 
than  other  methods. 

Consider  fitting  the  function 


Y  =  f(x;#)  (23) 

to  a  set  of  data  points  (x, ,  Y()  i  —  1,2, . .  .  L,  f(x;0)  nonlinear  in  the  parameters 
represented  by  the  vector  0=  (0o,#i . 02>»). 


Y,  -  f  (x,  ;*) 

Now,  using  ratios - to  meet  our  needs,  instead  of  the  customary 

f(xi  ;«) 

deviations  Y,  —  f(X(  ;<»),  we  want  to  minimize 


tr 


\ 


I 


L 

♦<•>  =  2  [<Y,  -  f(x,;*))/f<x,;#)]- 

i  i 

L 

=  S  tY,f-'(x,;#)  -  IP  . 

i  i 


Let 


K<x,;»)  =  . 


» 


(24) 


(25) 
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The  Taylor-series  expansion  of  g(xi  ;<?),  to  first  order  terms,  about  the  vector 
$w  =  (0<«>,  ew, . . .  ,0<°>),  vector  of  preliminary  parameter  estimates,  is 

0  1  2M 


2M 

g(x,;«)  =  g(x,;e<®>)  +  2  ) 

l=o 


(26) 


where 


Di  =  «i  -  «j‘0)  , 

dg(x,;«) 

g.(x,;»<®n  = - |  »  =  «<®> 

d«, 


Tliis 


Ll  £al  x 

*<*)  =  2  Y,[g(x, ;»(«>)  +  2  Djg.Cx,;^®))]  -  1  *  . 

1=1  l  l=o  ) 

Partial  derivative  of  <&(0)  with  respect  to  is 


(27) 


L  2M 

*,,<«)  =22  Y|[g(x1;*<®>)  +  2  DjgJ(X|;#<°))]  —  1  Y,gh  (x, ;#«») 

1=5  1=0 


(28) 


Setting  <M9)  =  0,  simplifying  and  transposing  terms,  we  have  the  normal 
equations 


2M 

2 

1=0 


2  (Y,gh  (x,;#<®>)  g^x,; 


1, 


,,<°>7]  Dj 


=  -  *  [Y,g(x,  ;«<»>)  -  l]Y,gh(x,;#<°>)  , 


(29) 


h  =  0,1 . 2M  , 

or  in  matrix  notation 


CD  =  E  . 


(80) 


Matrix  equation  80  is  solved  for  the  correction  vector,  D  which  in  turn 
is  used  in  obtaining  the  parameter  estimates. 


The  next  trial  vector  in  the  iterative  process  is 
=  »,<«>  +  Dj,  j  s  0,1, ...... 


(81) 


This  is  considered  the  Gauss  or  Gauss-Newton  approach  to  the  estima¬ 
tion  problem.  In  practice,  instead  of  using  Dj  as  is,  a  step  sise  FDj,  0  <  P  ^  1, 
is  used  in  an  effort  to  increase  convergence.  This  is  in  contrast  to  the  gradient 
methods  that  use  full  steps  in  the  direction  of  the  negative  gradient. 
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Proceeding  further  to  the  basic  construction  of  Marquardt’s  algorithm, 
equation  30  is  modified  by  adding  an  arbitrary  constant  A  to  each  diagonal  ele¬ 
ment  of  C.  Then  we  have  for  the  rth  iteration,  the  matrix  equation 

[C*(f>  +  \  <f)  I]  D*(0  =  E*(o  ,  (32) 

which  is  solved  for  the  column  matrix  D*<r),  where 

C*('>  =  [<:♦<')] 

uj 

=  [c(')/(c('))»(c(»)%]  ,  (33) 

h  j  hh  j  J 

E*(*)  =  [e*(r>] 

=  [ejr>  /  <•«)*]  ,  (34) 

and  chJ  and  e}  are  elements  of  C  and  E  of  equation  30. 

From  D*(r)  we  obtain 


D<r>  =  D*(f)  /  (c<'>)%  ,  j  =  0,1 . 2M  .  (35) 

i  1  )i 

The  next  vector  used  in  the  process  has  components 

#('+»)  =  *<0  +  D<r>  ,  j  =  0,1 . 2M  (36) 

The  choice  of  the  constant  A<r)  is  critical.  By  trial  and  error  it  is  deter¬ 
mined  so  that  $<r>  <  4><r  l». 


The  strategy  employed,  along  lines  somewhat  similar  to  those  laid  down 
by  Marquardt,  is  as  follows: 

Let  A(r  l)  denote  the  value  of  A  that  is  associated  with  the  (r-l)st  itera¬ 
tion,  initial  values  A(0>  =  10'*  and  =  (PMSR)(DF).  Choose  a  constant 
v  >  1,  a  suitable  value  chosen  is  10,  and 

Test  1.  If  A<,-,)  10_T,  proceed  to  test  2.  If  not,  let  A,r)  =  A,r-,>/»> 

and  compute  ♦<r)(A<,‘)- 

If  then  the  parameter  estimates  are  accepted; 

otherwise,  proceed  to  test  2. 

Teat  2.  Let  a*'*  =  A(,_l>  and  compute  ♦<r>(A<r)) . 


If  ♦<'>  <  then  the  parameter  estimates  are  accepted;  con- 
trarily,  set  m  —  1  and  proceed  to  test  3. 

Test  3.  If  a(’>  —  10l#,  the  process  is  considered  to  be  diverging  and 
further  cycling  is  terminated.  When  the  inequality  does  not  hold,  we  let 
a<m  ==  and  compute  ♦,r,(A'r>). 


If  «<r)  <  then  the  parameter  estimates  are  accepted; 

otherwise,  proceed  to  test  4.  It  should  be  observed  that  the  values  of  x<r>  can 
become  extremely  large  in  the  case  of  parameter  estimates  with  high  (>  .99) 
correlation. 

Test  4.  Compute  an  angle 

[2M  2M  2M  1 

2  d<r>e<r>  /  (2  d2<*>)K  (2  .  (37) 

J-0  J  i-0  1  J-0  1  _J 

If  y(f)  ^  y0  =  »/ 4,  then  increment  w  by  1  and  return  to  test  3. 

If  y<r)  <  y0  ,  then  proceed  to  find  a  constant  F<r>  such  that 
F(r>D(r>,  0  <  F<r)  ^  1,  is  the  step  vector;  i.e.,  we  want 

tW  =  e  (r-n  +  Fi'iD(f)  .  (38) 

This  constant  F(r)  is  found  by  raising  the  fraction  y%  to  succes¬ 
sive  powers  until  either  F(r)  ^  10-10,  in  which  case  «>  is  incremented  by  1  and  a 
return  to  test  3  is  made  after  setting  F(r>  =■■  1,  or  it  is  found  that  *(r)  < 
and  the  parameter  estimates  are  accepted. 

For  use  further  on  in  the  sequel,  we  have 
M 

g(X,;«)  =  [<*0  +  2  On.  «*P<  -  . 


=  s-*  (X,;a^)  , 

(39) 

g„  (x,;«)  =  -  S--‘  (Xpa./OS,,  (x,;aji)  , 

(40) 

h  =  0,1 . 2M  , 

where 

/I  ,h=0 

S„  (X,;«jj)  =  |  “p<  “  ^o+iX|)  ’  h  =  1,S*  • 

. . ,  odd  integer 

\  —  •k/gXftxpi  —  p  1,/jX,)  ,  h  =  2,4, . 

. . , even  integer 

and  a,  p  are  vectors  defined  as : 

•  —  (*0>  •» . . )  • 

(•  —  (J|.  4|>  ■  •  •  >0||)  • 


Hence  the  normal  equations  29  become 


SM 

X 

i-o 


L.  . 


t-t 


Dj 


L 

=  X  [Y,8-Mxt;«'“»jf"»)  -  l)Y,8-s(xl;«<#»rfM**)Stl(xl;««*>d»‘B»)  . 
h  =  0.1 . SH  . 


(41) 


Computer  version 

We  continue  the  sequel  with  a  description  of  the  computer  program  based  on 
Marquardt’*'  algorithm.  This  description  will  include  program  input,  general 
computation?’  steps,  additional  computations,  and  output.  Then  follows,  for 
completeness,  a  program  flow  chart,  description  of  flow  chart,  subprogram  usage, 
computer  memory  requirements,  and  several  numeric  examples  of  analyses  of 
both  simulated  and  empirical  data. 


Input: 

1.  '  Program  parameters. 

a.  From  a  punched  card. 

NIA:  Maximum  number  of  iterations  allowed. 

b.  From  output  tape  produced  by  “peel-off”  program. 

(1)  EXP:  Experiment  number. 

(2)  SAMP:  Sample  number. 

(3)  MODEL:  1  for  no  asymptote, 

2  otherwise. 

(4)  L:  Sample  size. 

(5)  M:  Number  of  components  to  refine. 

(6)  NP :  Number  of  parameters  to  refine. 

(7)  PSMR:  “Unrefined"  mean  square  ratio  based  on  preliminary 
parameter  estimates. 

2.  Preliminary  estimates  of  the  model  parameters  (from  output  tape). 

»  =  U . M. 


3.  Data  (from  output  tape), 
x,  Y  values. 

General  computational  steps: 

Step  1.  Compute  (initial)  gage  quantity 

=  (PM8RHDF)  , 


where 


DF  =  L  -  811  for  Modol  1 

=  L  —  (2M  +  1)  tor  Mofel  8. 


(48) 
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Step  2.  Compute  the  elements  of  the  matrices  C*<'>,  E*<r>, 

%  V4 

C*<r>  =  C<r>  /  (C<r»)  (c(f))  ,  (48) 

hj  hi  hh  Ji  '  ; 

where 

L  , 

c<r>  =  2  YtS-»(x1;«<r>^<'>)Sil(x1;«<')^f))S)(*l;«<f»^))  , 

hj  1-1 

h,j  =  0,1.- ...2M  , 

and  S(x,  \a(r) ,p(T)) ,  Sh(xi  ;a(r),/8<r>)  are  as  defined  in  equations  39,  40. 

e*(r)  =  e<r>  /  <c('))W  ,  (44) 

i  j  ii 

where 

L 

e<r>  =  2  [Y,S-Mx,;a<r> ,£<'>)  _  l]Y,S-2(xi;«('>^('))  SAx^'K^)  , 

i  i-i 

j  =  0,1, ...,2M  . 

Step  3.  Solve  matrix  equation 

[C*<r>  +  \<OI]D*<r>  =  E*(')  .  (45) 

Step  4.  Apply  strategy,  outlined  by  series  of  tests  given  above,  for  choosing 
A.(r)  and  selecting  F<r). 

Step  5.  Terminate  process  in  accordance  with  tests  : 

|Dj<r>|/(r  +  |»,<») |)  <  «  ,  j  =  0 . 2M  (46) 

where 

»(r>  =  <«(*>  ,  «<r>  ,  «<'>  ,  «<n  .  0<J>)  , 

r  ;  arbitrary,  say  10  ~s  , 

«  arbitrary,  say  5(10~T)  , 


or 


N  =  NIA 


(47) 


where 


N  Number  of  completed  iterations, 

NIA:  Maximum  number  of  iterations  allowed. 

Step  6.  Return  to  step  2,  if  neither  test  in  step  5  holds,  with  the  current 
solution  vector 

#(M  =  <»<«  ,  *<»>  ,  . •<'»  ,  0<J>)  (48) 

as  the  initial  vector  along  with  A<r>  value  for  the  next  iteration. 


25 


Additional  computations: 

1.  Compute  omegas 

«„,(»>  =  exp(  -  pmM)  ,  m  =  1,2,...,M  , 

where 

0m<B>  =  Final  value  of  p  for  mth  component. 

2.  Compute  estimated  Y  values 

M 

YEST.  =  «0(»)  +  2  am<“>  exp(  -  0m<n>x,>  ,  i  =  1,2, . . . ,  L  . 

m— 1 

3.  Compute  ratios 

RA,  =  (Y,  -  YEST,)  /  YEST,  ,  1  =  1*..  .L. 


4.  Compute  cumulative  sums  of  squared  ratios 


CRATIO,  =  2  RAf  ,  i  =  1* . L  . 

J-i 


5.  Compute  “refined”  mean  square  ratio 
FMSR  =  CRATIO,/DF  , 

where 

DF  =  L-2M  for  Model  1 

=  L  —  (2M  +  1)  for  Model  2  . 


Output: 

1.  Experiment  number. 

2.  Sample  number 
S.  Data. 

4.  Estimated  Y  values. 


(49) 


(50) 


(61) 


(52) 


(53) 
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6.  Ratios. 

6.  Cumulative  sums  of  squared  ratios. 


7.  Model. 


8.  Number  of  data  points. 

9.  Number  of  iterations  taken. 

10.  Number  of  iterations  allowed. 

11.  "Unrefined”  mean  square  ratio. 

12.  “Refined”  mean  square  ratio. 

13.  Preliminary  and  final  estimates  of  constant,  alphas,  betas,  and  omegas. 
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FLOW  CHART  I 

FLOW  CHART  OF  THE  COMPOSITE  GAUSS- NEWTON  AND  GRADIENT  METHOD 


FLOW  CHART  K  CONT. 


l**>.llr,*i  o.  **|  J 

WAi  «l. HP  1  1 


it 


«CT  ITDIT 


Description  of  flow  chart  II 


Box  1  -  S:  Self-explanatory. 

4:  The  following  arrays  are  cleared — 
a:  Parameter  array. 
p:  Parameter  array. 
u :  Parameter  array. 

(n— 1) at  parameter  iterant  array. 
e <»> :  nth  parameter  iterant  array. 

D:  Parameter  step  adjustments  matrix. 

E:  A  basic  column  matrix  used  in  the  iterative  process.  Its  companion  matrix 
is  basic  matrix  C. 

S:  Square  roots  of  the  diagonal  elements  of  matrix  C. 

ES:  Column  matrix  produced  by  dividing  elements  of  matrix  E  by  the  cor¬ 
responding  elements  of  S-array. 

T:  First  partial  derivatives  of  the  mathematical  model  with  respect  to  its 
parameters. 

X:  Abscissae  of  the  data  points. 

Y :  Ordinates  of  the  data  points. 

5  -  7 :  Self-explanatory. 

8:  Test  parameter  X  .  If  the  inequality  holds,  go  to  TEST  2 — box  12. 

9:  Self-explanatory. 

10:  See  description  of  subroutine  GET  ♦,  boxes  23-38. 

11:  Self-explanatory. 

12:  See  description  of  subroutine  GET  4,  boxes  23-38. 

18:  Test  parameter  X  .  If  the  inequality  holds,  go  to  DIVERGE — box  22. 

14:  Self-explanatory. 

18:  See  description  of  subroutine  GET  0,  boxes  23-38. 

10-18:  Compute  an  angle  y  .  If  the  inequality  holds,  go  to  TEST  8— box  18. 

19:  Self-explanatory. 

29:  Teat  step  adjustment  factor  F.  If  the  Inequality  holds,  go  to  TEST  8 — box  18. 

21 :  See  description  of  subroutine  TEST  4,  boxes  89-46. 

22:  Self-explanatory. 

28:  Self-explanatory. 

24:  Clear  locations  for  elements  of  matrices  C,  E,  ES,  S  and  T. 

28:  Compute  first  derivative  of  tbs  function  in  the  model  with  respect  to  the 
constant  tarn.  Counter  k  la  aet  equal  to  one. 

28-28:  Construction  of  basic  matrices  C  and  E. 

28:  Self  -explanatory . 

28:  Teat  counter  k  against  gage  word  L  If  the  inequality  holds,  go  to  GET  C  and  E — 

box  26. 
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31 :  Compute  adjustment  quantities  for  matrices  C  and  E. 

32:  Transform  matrix  C  into  simple  correlation  matrix,  resulting  matrix  again 
called  C;  adjust  matrix  E,  resulting  matrix  called  ES. 

33:  Add  X  to  the  main  diagonal  elements  of  simple  correlation  matrix  C. 

34:  Compute  inverse  of  modified  simple  correlation  matrix  C. 

35:  Compute  product  matrix  D,  D  =  C-1  ES. 

36:  Transform  elements  of  matrix  D*  bank  into-original  units  of  measurement. 

37 :  fHf-explanatory. 

38:  See  description  of  subroutine  TEST  *,  boxes  39-46. 

39:  Self-explanatory. 

40:  Compute  the  nth  iterant  of  parameter  estimates. 

41:  Evaluate  nth  iterant  *  function. 

42:  Test  (n— l)st  and  nth  iterant  ♦  functions.  If  the  inequality  holds,  return  to 
location  specified  by  subroutine  exit. 

43:  Test  parameter  estimates  for  acceptability.  If  the  inequality  fails  to  hold,  go  to 
OUTPUT— box  47. 

44:  Self-explanatory. 

45:  Test  number  of  iterations  completed  counter  n.  If  the  inequality  does  not  hold, 
go  to  OUTPUT— box  47. 

46:  Initialize  »("-i)-array  for  next  iteration,  go  to  TEST  1 — box  8. 

47:  Relocation  of  the  final  estimate  of  «o»  a>s»  and  0’s;  computation  of  the  associated 
w  estimates. 

48-49:  Print  output: 

a.  EXP.  SAMP. 

b.  X,  Y,  YE,  RA,  CRATIO  for  each  data  point. 

c.  EXP,  SAMP,  MODEL,  L,  N,  NIA,  PMSR,  PMSR. 

d.  a0<°>  and  estimates  «<°>,  0<">,  «<°>  for  each  component  (preliminary 
estimates). 

e.  •**<">  and  estimates  «(»),  0<«i,  „<«>  for  each  component  (final  estimates). 

Subprogram  usage 

Function  subprograms  (the  first  one  is  a  standard  library  function)  listed 
below  proved  very  helpful: 

1.  EXPF.  Argument:  A  (location  of  expression  A).  Function:  computes 
the  value  exp  (A). 

2.  YEST.  Argument:  K  (location  of  subscript  K).  Function:  computes 
the  kth  estimated  Y  value  using  the  parameter  estimates  a*">  ,  a"*  ,  and  £<*>  , 

0  I  I 

i  —  1,2, . . .  ,M;  n  denotes  nth  iteration. 

Explicitly  presented  in  the  flow  chart  are  two  subroutine  subprograms: 

1.  GET  ♦.  Function :  computes  full  step  adjustments  for  the  parameter 
estimates  under  refinement. 
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2.  TEST  $.  Function:  adjusts  the  parameter  estimates  using  a  full  or  frac¬ 
tional  part  of  the  step  adjustments  and  then  computes  the  sum  of  squared 
ratios.  Performs  a  number  of  tests:  *<">  tested  against  determines  if 

parameter  estimates  produced  during  the  nth  iteration  meet  the  acceptance  crite¬ 
rion,  or  if  permissible  number  of  iterations  criterion  is  met  in  case  of  failure 
of  parameter  estimates  to  satisfy  acceptance  criterion. 

Two  standard  library  subroutine  subprograms  not  explicitly  shown  in  the  flow 
chart  but  of  paramount  importance  in  the  iterative  process  are: 

1.  AMMAIN.  Arguments:  DET  (location  for  value  of  the  determinant  of 
the  matrix),  N  (location  for  dimension  of  the  matrix),  C  (location  of  the  first 
word  of  the  matrix  to  be  inverted),  and  T  (location  of  the  first  of  2N  words  used 
as  temporary  storage).  Function:  computes  the  inverse  and  determinant  of  a 
matrix  C  in  single-precision  floating-point  arithmetic.  The  Gauss  process  of 
elimination  is  used. 


2.  AMMAMU.  Arguments:  NR  (location  for  row  dimension  of  matrix  C), 
NC  (location  for  column  dimension  of  matrix  C),R  (location  for  column  dimen¬ 
sion  of  matrix  E),  G  (location  of  first  element  of  matrix  C),  E  (location  of  first 
element  of  matrix  E),  D  (location  of  first  element  of  product  matrix  D).  Func¬ 
tion  :  forms  the  matrix  product  D  =  CE. 


Memory  requirements 


Program 

About  3,080  words 

a<°) 

0 

1  word 

«(«) 

i 

M  words 

»<«> 

i 

M  words 

*,«») 

i 

M  words 

0<—i) 

t 

2M  +  1  words 

i 

2M  f  1  words 

„<»> 

i 

M  words 

X -array 

L  words 

Y-array 

L  words 

C-array 

(2M  + 1)*  words 

D-array 

2M  + 1  words 

E-array 

2M  4- 1  words 

ES-array 

2M  4  1  words 

Other 

170  words 

Total:  3,256  +  14M  +  2L  +  (2M  4- 1)*  words  approxima 
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Examples 

Parameter-estimate  input  to  the  composite  Gauss-Newton  and  gradient  meth¬ 
od  computer  program  consisted  of  the  preliminary  estimates  shown  in  tables 
I  to  VI.  The  refined  parameter  estimates  along  with  supplementary  data  are 
presented  in  tables  VII  to  XII.  Here  results  for  simulated  data  are  reserved  for 
tables  VII  to  IX  while  for  empirical  data  tables  X  to  XII  are  utilized. 


TABLE  VII 

Refinement  of  parameter  estimates 
(Simulated  data) 

Sample  size:  50 
Value  of  p:  0.0005 
Number  of  iterations:  5 

Square  root  of  MS  ratio  for  “unrefined”  parameter  estimates:  0.00725 
Square  root  of  MS  ratio  for  “refined”  parameter  estimates:  0.00047 


Parameter 

True  value 

Preliminary 

estimate 

Final  estimate 

•o 

0.00500 

0.00600* 

0.00497 

0.02000 

0.02047 

0.02000 

Pi 

0.03046 

0.08561 

0.03034 

“1 

0.97000 

0.96602 

0.97012 

■* 

0.70000 

0.69812 

0.69981 

Pt 

0.35667 

0.85739 

0.35669 

“a 

0.70000 

0.69950 

0.70006 

i 

'BrtimaU  of  a,  arbitrarily  choaaa  in  tha  naighborbood  of  tbo  true  nht. 
Running  tlman— 


Preliminary :  0.1  min. 

Final  1.7  min. 


The  data  points  listed  below  were  used  to  derive  the  results  presented 
in  table  VII;  thus,  they  may  serve  as  test  data. 


Test  Data 


X 

1 

2 

8 

4 

6 

Y 

0.614116904 

0J66671076 

0.243296656 

0.190018417 

0.189906807 

l 

6 

7 

8 

9 

10 

Y 

0.104017778 

0.078778202 

0.061029668 

0.048484872 

0.089610206 

X 

11 

12 

18 

14 

16 

Y 

0.033133061 

0.028666668 

0.026232096 

0.022818472 

0.020982678 

X 

16 

17 

18 

19 

20 

Y 

0.019616327 

0.018663200 

0.017690670 

0.016996348 

0.016446828 

X 

21 

22 

23 

24 

26 

Y 

0.016943802 

0.016604852 

0.016120464 

0.014761646 

0.014443461 

X 

26 

27 

28 

29 

30 

Y 

0.014126197 

0.013832992 

0.013563697 

0.013294183 

0.013084266 

X 

31 

32 

33 

34 

36 

Y 

0.012786063 

0.012650016 

0.012322636 

0.012106027 

0.011888306 

X 

36 

37 

38 

39 

40 

Y 

0.011691969 

0.011487734 

0.011281136 

0.011103610 

0.010918208 

X 

41 

42 

43 

44 

46 

Y 

0.010739666 

0.010669082 

0.010402436 

0.010236416 

0.010073497 

X 

46 

47 

48 

49 

60 

Y 

0.009914986 

0.009782278 

0.009636448 

0.009496308 

0.009359069 

TABLE  VIII 

Refinement  of  parameter  estimates 
(Simulated  data) 

Sample  aiie:  76 
Value  of  p:  0.0010 
Number  of  iterations:  6 

Square  root  of  MS  ratio  for  “unrefined”  parameter  estimates:  0.0433 
Square  root  of  MS  ratio  for  “refined”  parameter  estimates:  0.0008 


Parameter 

True  value 

Preliminary 

estimate 

Final  estimate 

•o 

None 

None 

None 

•i 

0.0400 

0.0649 

0.0400 

fix 

0.0806 

0.0868 

0.0806 

“I 

0.9700 

0.9648 

0.9700 

•i 

0.4000 

0.8977 

0.3997 

fit 

0.0943 

0.1029 

0.0948 

0.9100 

0.9022 

0.9100 

•* 

0.2000 

0.1796 

0.1987 

fit 

0.4806 

0.4812 

0.4270 

•» 

0.6600 

0.6180 

0.6625 

Resales  t law 

FnUmlaaryi  M  at*. 
riMl  t  L4  ale. 


i 

i 
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TABLE  IX 


Refinement  of  parameter  estimates 
(Simulated  data) 

Sample  size:  42 
Value  of  p:  0.0030 
Number  of  iterations:  89 

Square  root  of  MS  ratio  for  "unrefined*'  parameter  estimates:  0.0125 
Square  root  of  MS  ratio  for  “refined”  parameter  estimates:  0.0028 


Parameter 

True  value 

Preliminary 

estimate 

Final  estimate 

«# 

0.0500 

0.0500* 

0.0500 

“1 

0.0500 

0.0243 

0.0474 

Pi 

0.1054 

0.0845 

0.1028 

»i 

0.9000 

0.9190 

0.9023 

“a 

0.1000 

0.0419 

0.0549 

Pi 

0.3667 

0.1622 

0.2976 

"a 

0.7000 

0.8503 

0.7426 

°s 

0.2000 

0.0061 

0.2121 

Ps 

0.6932 

0.2055 

0.5723 

"3 

0.5000 

0.8143 

0.5642 

°4 

0.4000 

0.5314 

0.4440 

P 4 

1.2040 

0.7203 

1.2233 

"4 

0.3000 

0.4866 

0.2943 

"Estimate  of  a0  simply  chosen  equal  to  the  true  value. 
Running  times — 

Preliminary:  1.2  min. 

Final  :  10.5  min. 


TABLE  X 

Refinement  of  parameter  estimates 
(Empirical  data) 

Sample  size:  29 
Number  of  iterations:  82 

Square  root  of  MS  ratio  for  "unrefined”  parameter 
estimates:  0.0876 

Square  root  of  MS  ratio  for  "refined”  parameter 
estimates:  0.0278 


Parameter 

Preliminary 

estimate 

Final  estimate 

•o 

0.0050* 

0.0062 

•i 

0.0119 

0.5566 

Pi 

0.0797 

0.2838 

-i 

0.9234 

0.7680 

«* 

0.8899 

0.3864 

P» 

0.3582 

0.5698 

**» 

0.6989 

0.5657 

of  a,  buni  on  t  vUunl  lupeetlon  of  pkHUd  eurvo. 
Running  Umea — 


Preliminary:  e.l  min. 
Klnai  3.4  min. 


TABLE  XI 


Refinement  of  parameter  estimates 
( Empirical  data) 

Sample  size:  49 
Number  of  iterations:  11 

Square  root  of  MS  ratio  for  “unrefined”  parameter 
estimates:  0.0665 

Square  root  of  MS  ratio  for  “refined”  parameter 
estimates:  0.0218 


Parameter 

Preliminary 

estimate 

Final  estimate 

«.> 

None 

None 

ai 

0.0147 

0.0107 

Pi 

0.0139 

0.0085 

“i 

0.9862 

0.9916 

a2 

0.1717 

0.1037 

Pi 

0.1625 

0.1143 

Wo 

0.8500 

0.8920 

a3 

0.5173 

0.6490 

Pi 

0.3623 

0.3922 

"3 

0.6961 

0.6766 

Running  times — 

Preliminary:  0.1  min. 
Final  :  1.5  min. 


TABLE  XII 

Refinement  of  parameter  estimates 
(Empincal  data) 

Sample  size:  140 
Number  of  iterations:  15 

Square  root  of  MS  ratio  for  “unrefined”  parameter 
estimates:  0.5504 

Square  root  of  MS  ratio  for  "refined”  parameter 
estimates:  0.C216 


Parameter 

Preliminary 

estimate 

Final  estimate 

•o 

0.0100* 

0.0116 

•» 

0.0047 

0.0908 

Pi 

0.0049 

0.0410 

"t 

0.9951 

0.9699 

•» 

0.0132 

0.1695 

Pi 

0.0279 

0.1477 

"i 

0.9725 

0.8627 

•» 

0.6208 

0.4374 

Pi 

0.2600 

0.4729 

0.7710 

0.6282 

‘btimate  of  a,  bated  on  a  vltual  latpoctioa  of  plotted  curv*. 
Rossini  ttmte-  ■ 

ProUBBlaory:  0.1  mla. 
rtnol  >  44  min. 


V.  COMMENTS 


Even  a  cursory  examination  of  the  basic  steps  pertinent  to  the  “peel-off* 
method  reveals  that  these  steps  are  quite  simple  in  mathematical  content  and  do 
not  offer  any  particular  problem  in  programming.  The  overall  “sticky”  areas 
of  “when  or  where  to  continue  or  discontinue”  the  fitting  process  no  doubt  can 
be  improved  upon.  What  and  how  much  improvement  is  a  question  that  can 
best  be  answered  through  extensive  testing  of  the  program  on  data  frequently 
handled  at  a  particular  computing  installation  or  laboratory.  Needless  to  say,  one 
should  obtain  the  best  possible  preliminary  estimate  of  the  constant  term  when 
such  is  present  in  the  model  being  fitted.  A  poor  preliminary  estimate  of  the 
constant  term  results  in  inferior  component  estimates,  and  this  in  turn  can 
result  in  a  great  increase  in  the  number  of  iterations  for  refinement  or  the  diver¬ 
gence  of  the  refining  process.  A  number  of  methods  for  constant-term  estima¬ 
tion  by  computer  were  tried  but  results  were  in  general  disappointing.  Further 
work  along  these  lines  is  contemplated.  Also  some  Monte  Carlo  studies  of  the 
parameters  estimates  are  planned.  All  raw  data  that  is  not  monotone  decreas¬ 
ing  should  be  smoothed  before  being  analyzed  and  preliminary  parameters  esti¬ 
mates  should  always  be  refined.  Lastly  the  use  of  the  upper  1%  points  of  the 
Snedecor-Fisher  (F)  distribution  in  the  F  test  as  proposed  and  used  by  the 
authors  may  provide  a  test  that  is  too  stringent  for  some  users.  In  such  cases 
the  appropriate  percentage  points  can  be  selected  on  the  basis  of  program  yield 
on  data  commonly  handled. 

One  should  be  able  to  program  the  composite  Gauss-Newton  and  gradient 
iterative  method  for  nonlinear  parameter  estimation  with  relative  ease.  The 
method  is  straightforward,  requiring  no  “special”  coding  technics.  Minor  changes 
envisaged  in  a  few  instances  are  those  relative  to  the  parameter  values  of  e, 
yo,  A(0),  »•,  and  r.  These  changes  may  be  made  in  accordance  with  the  immediate 
requirements  of  the  user.  In  passing,  one  might  consider  the  use  of  double-preci¬ 
sion  floating  point  in  computing  *  with  Y,  in  single  or  double  precision  for 
tests  1  and  2  (see  strategy  for  choosing  A)  if  round-off  contributes  to  erratic 
fluctuations  in  <t> . 
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